#ifndef FF_NN2_SSE_COR_TH_N_H
#define FF_NN2_SSE_COR_TH_N_H

#include <math.h>
#include <xmmintrin.h>
#include <iostream>


// my includes
#include "utils/utils.02.h"


namespace dj01 {


class FJPF {
	
	public:
	
	uint n, N;
	double *m;
	double *x0, *x1, *x2;
	double *v0, *v1, *v2;
	double *f0, *f1, *f2;
	double *j0, *j1, *j2;
	double *pot;
	double eps2;
	uint n_core;
	
	lluint nij, nij_max, nij_min; // work length for [1, nth)-threads and nth-thread
	Timer tr;
	
	// n, j1 should be even, j2 should be odd, arra
	inline virtual void ff(	uint n, uint N, 
									uint i1_,uint j1_,uint i2_,uint j2_, 
									double*m, 
									double*x0,double*x1,double*x2,
									double*v0,double*v1,double*v2,  
									double*f0,double*f1,double*f2, 
									double*j0,double*j1,double*j2, 
									double*pot, double eps2) {

		if (n%2!=0 && n>1) n++; // we need even number for SSE
		
		uint dj=10000; // array indexing range (must be align 2) 10000 is optimal
		uint kn = ceil((double)N/dj);
		uint _i1, _i2, _j1, _j2, _j11, _j21, _j12, _j22, jL, jR ;	
	
		double tmp[2] __attribute__((aligned(16)));	

		__m128d _eps2 =_mm_setr_pd(eps2, eps2);
		__m128d _3 =_mm_setr_pd(3, 3);
		
		
		uint _N=N+1;
		
		//memset(m, 0, (N+1)*sizeof(*m));
		
		#ifndef ACC_OFF
		memset(f0, 0, _N*sizeof(*f0));
		memset(f1, 0, _N*sizeof(*f1));
		memset(f2, 0, _N*sizeof(*f2));
		#endif
		
		#ifndef JERK_OFF
		memset(j0, 0, _N*sizeof(*j0));
		memset(j1, 0, _N*sizeof(*j1));
		memset(j2, 0, _N*sizeof(*j2));
		#endif
		
		#ifndef POT_OFF
		memset(pot, 0, _N*sizeof(*pot));
		#endif
		
		for (uint k=0; k<kn; k++) {  // lets devide j walk on kn parts with size dj performance reasons

			jL= k*dj; // lets define Light j-coordinate 
			N-jL>dj ? jR=jL+dj-1 : jR = N-1; // j-Left coordinate
				
			_i1=i1_; _i2=i2_; // ??? should be modified later on  ???
			
			for (uint i=_i1; i<=_i2; i++) {
				
				__m128d mi =_mm_setr_pd(m[i], m[i]);
				
				__m128d x0i=_mm_setr_pd(x0[i], x0[i]);
				__m128d x1i=_mm_setr_pd(x1[i], x1[i]);
				__m128d x2i=_mm_setr_pd(x2[i], x2[i]);
				
				#ifndef ACC_OFF
				__m128d f0i=_mm_setr_pd(0, 0);
				__m128d f1i=_mm_setr_pd(0, 0);
				__m128d f2i=_mm_setr_pd(0, 0);
				#endif
				
				#ifndef JERK_OFF
				__m128d v0i=_mm_setr_pd(v0[i], v0[i]);
				__m128d v1i=_mm_setr_pd(v1[i], v1[i]);
				__m128d v2i=_mm_setr_pd(v2[i], v2[i]);
				
				__m128d j0i=_mm_setr_pd(0, 0);
				__m128d j1i=_mm_setr_pd(0, 0);
				__m128d j2i=_mm_setr_pd(0, 0);
				#endif
				
				#ifndef POT_OFF
				__m128d poti=_mm_setr_pd(0, 0);
				#endif
				
				i==i1_ ? _j1=j1_ : _j1=i+1;
				if ((_j1==i+1) && (i%2==0)) _j1=i;
				if (_j1<jL) _j1=jL;
		
				i==i2_ ? _j2=j2_ : _j2=N-1;
				if (jR<_j2) _j2=jR;
				
				if (_j1>_j2) continue;
				if (_j2%2==0) _j2++;
					
				// to spread j tasks among nn2 and nn part
				if (_j1<n-1) {
					
					_j11=_j1; 
					if (n-1<_j2){ _j12=n-1; _j21=n; _j22=_j2; } 
							 else	{ _j21=_j22=0; _j12=_j2;      }	
				
				}  else { 
					
					_j11=_j12=0;	
					_j21=_j1; _j22=_j2; 
					
				}
				
				// nn2
				for (uint j=_j11; j<_j12; j+=2) {
							
					__m128d x0j=_mm_load_pd(&x0[j]); 
					__m128d dx0ij=_mm_sub_pd(x0i, x0j);
					__m128d sqrdxij0=_mm_mul_pd(dx0ij, dx0ij);
					
					__m128d x1j=_mm_load_pd(&x1[j]);
					__m128d dx1ij=_mm_sub_pd(x1i, x1j);
					__m128d sqrdxij1=_mm_mul_pd(dx1ij, dx1ij);
					
					__m128d x2j=_mm_load_pd(&x2[j]);	
					__m128d dx2ij=_mm_sub_pd(x2i, x2j);
 					__m128d sqrdxij2=_mm_mul_pd(dx2ij, dx2ij);

					#ifndef JERK_OFF
 					__m128d v0j=_mm_load_pd(&v0[j]); 
 					__m128d dv0ij=_mm_sub_pd(v0i, v0j);
 					
					__m128d v1j=_mm_load_pd(&v1[j]);
					__m128d dv1ij=_mm_sub_pd(v1i, v1j);
					
					__m128d v2j=_mm_load_pd(&v2[j]);	
					__m128d dv2ij=_mm_sub_pd(v2i, v2j);
					#endif
					
					__m128d rij2 = _mm_add_pd(sqrdxij0, sqrdxij1);
					rij2 = _mm_add_pd(rij2, sqrdxij2);
					rij2 = _mm_add_pd(rij2, _eps2);
			
					__m128d rij = _mm_sqrt_pd(rij2);
					__m128d rij3 = _mm_mul_pd(rij2, rij);

					__m128d mj=_mm_load_pd(&m[j]);
					__m128d mimj = _mm_mul_pd(mi, mj);
					__m128d mimjrij3 = _mm_div_pd(mimj, rij3);
					
					#ifndef JERK_OFF
					__m128d dxijdvij0 =_mm_mul_pd(dx0ij, dv0ij);
					__m128d dxijdvij1 =_mm_mul_pd(dx1ij, dv1ij);
					__m128d dxijdvij2 =_mm_mul_pd(dx2ij, dv2ij);
					
					__m128d 
					xv3rij2 = _mm_add_pd(dxijdvij0, dxijdvij1);
					xv3rij2 = _mm_add_pd(xv3rij2, dxijdvij2);
					xv3rij2 = _mm_mul_pd(xv3rij2, _3);
					xv3rij2 = _mm_div_pd(xv3rij2, rij2);
					#endif
					
					// POT 
					#ifndef POT_OFF
					__m128d potij = _mm_mul_pd(mimjrij3, rij2);
 					poti =_mm_sub_pd(poti, potij);	
					__m128d potj =_mm_load_pd(&pot[j]);
					potj =_mm_sub_pd(potj, potij);
					_mm_store_pd(&pot[j], potj);
					#endif
					
					#ifndef ACC_OFF
					__m128d f0ij = _mm_mul_pd(mimjrij3, dx0ij);
					f0i =_mm_sub_pd(f0i, f0ij);
					__m128d f0j=_mm_load_pd(&f0[j]);
					f0j =_mm_add_pd(f0j, f0ij);
					_mm_store_pd(&f0[j], f0j);
					
					__m128d f1ij = _mm_mul_pd(mimjrij3, dx1ij);
					f1i =_mm_sub_pd(f1i, f1ij);
					__m128d f1j=_mm_load_pd(&f1[j]);
					f1j =_mm_add_pd(f1j, f1ij);
					_mm_store_pd(&f1[j], f1j);
					
					__m128d f2ij = _mm_mul_pd(mimjrij3, dx2ij);
					f2i =_mm_sub_pd(f2i, f2ij);
					__m128d f2j=_mm_load_pd(&f2[j]);
					f2j =_mm_add_pd(f2j, f2ij);
					_mm_store_pd(&f2[j], f2j);
					#endif
					
					#ifndef JERK_OFF
					__m128d tmp0 =	_mm_mul_pd(xv3rij2, dx0ij);
					tmp0 =	_mm_sub_pd(dv0ij, tmp0);
					__m128d j0ij = _mm_mul_pd(mimjrij3, tmp0);
					j0i =_mm_sub_pd(j0i, j0ij);
					__m128d j0j=_mm_load_pd(&j0[j]);
					j0j =_mm_add_pd(j0j, j0ij);
					_mm_store_pd(&j0[j], j0j);
					
					__m128d tmp1 =	_mm_mul_pd(xv3rij2, dx1ij);
					tmp1 = _mm_sub_pd(dv1ij, tmp1);
					__m128d j1ij = _mm_mul_pd(mimjrij3, tmp1);
					j1i =_mm_sub_pd(j1i, j1ij);
					__m128d j1j=_mm_load_pd(&j1[j]);
					j1j =_mm_add_pd(j1j, j1ij);
					_mm_store_pd(&j1[j], j1j);
					
					__m128d tmp2 =	_mm_mul_pd(xv3rij2, dx2ij);
					tmp2 = _mm_sub_pd(dv2ij, tmp2);
					__m128d j2ij = _mm_mul_pd(mimjrij3, tmp2);
					j2i =_mm_sub_pd(j2i, j2ij);
					__m128d j2j=_mm_load_pd(&j2[j]);
					j2j =_mm_add_pd(j2j, j2ij);
					_mm_store_pd(&j2[j], j2j);
					#endif
 					
				}

				#ifdef DEBUG
				printf("j2=[%u, %u]\n", _j21, _j22);
				#endif
				
				//nn
				for (uint j=_j21; j<_j22; j+=2) {
					
					__m128d x0j=_mm_load_pd(&x0[j]); 
					__m128d dx0ij=_mm_sub_pd(x0i, x0j);
					__m128d sqrdxij0=_mm_mul_pd(dx0ij, dx0ij);
 										
					__m128d x1j=_mm_load_pd(&x1[j]);
					__m128d dx1ij=_mm_sub_pd(x1i, x1j);
					__m128d sqrdxij1=_mm_mul_pd(dx1ij, dx1ij);

					__m128d x2j=_mm_load_pd(&x2[j]);	
					__m128d dx2ij=_mm_sub_pd(x2i, x2j);
 					__m128d sqrdxij2=_mm_mul_pd(dx2ij, dx2ij);
					
					#ifndef JERK_OFF
					__m128d v0j=_mm_load_pd(&v0[j]); 
 					__m128d dv0ij=_mm_sub_pd(v0i, v0j);
					
					__m128d v1j=_mm_load_pd(&v1[j]);
					__m128d dv1ij=_mm_sub_pd(v1i, v1j);
					
					__m128d v2j=_mm_load_pd(&v2[j]);	
					__m128d dv2ij=_mm_sub_pd(v2i, v2j);
					#endif
			
					__m128d rij2 = _mm_add_pd(sqrdxij0, sqrdxij1);
					rij2 = _mm_add_pd(rij2, sqrdxij2);
					rij2 = _mm_add_pd(rij2, _eps2);
			
					__m128d rij = _mm_sqrt_pd(rij2);
					__m128d rij3 = _mm_mul_pd(rij2, rij);

					__m128d mj=_mm_load_pd(&m[j]);
					__m128d mimj = _mm_mul_pd(mi, mj);
					__m128d mimjrij3 = _mm_div_pd(mimj, rij3);
					
					// POT 
					#ifndef POT_OFF
					__m128d potij = _mm_mul_pd(mimjrij3, rij2);
					poti =_mm_sub_pd(poti, potij);	
					#endif

					
					#ifndef ACC_OFF
					__m128d f0ij = _mm_mul_pd(mimjrij3, dx0ij);
					f0i =_mm_sub_pd(f0i, f0ij);
					
					__m128d f1ij = _mm_mul_pd(mimjrij3, dx1ij);
					f1i =_mm_sub_pd(f1i, f1ij);
										
					__m128d f2ij = _mm_mul_pd(mimjrij3, dx2ij);
					f2i =_mm_sub_pd(f2i, f2ij);
					#endif					

					#ifndef JERK_OFF
					__m128d dxijdvij0 =_mm_mul_pd(dx0ij, dv0ij);
					__m128d dxijdvij1 =_mm_mul_pd(dx1ij, dv1ij);
					__m128d dxijdvij2 =_mm_mul_pd(dx2ij, dv2ij);
					
					__m128d 
					xv3rij2 = _mm_add_pd(dxijdvij0, dxijdvij1);
					xv3rij2 = _mm_add_pd(xv3rij2, dxijdvij2);
					xv3rij2 = _mm_mul_pd(xv3rij2, _3);
					xv3rij2 = _mm_div_pd(xv3rij2, rij2);
					
					__m128d 
					tmp0 = _mm_mul_pd(xv3rij2, dx0ij);
					tmp0 = _mm_sub_pd(dv0ij, tmp0);
					__m128d j0ij = _mm_mul_pd(mimjrij3, tmp0);
					j0i =_mm_sub_pd(j0i, j0ij);
					
					__m128d 
					tmp1 = _mm_mul_pd(xv3rij2, dx1ij);
					tmp1 = _mm_sub_pd(dv1ij, tmp1);
					__m128d j1ij = _mm_mul_pd(mimjrij3, tmp1);
					j1i =_mm_sub_pd(j1i, j1ij);
					
					__m128d 
					tmp2 = _mm_mul_pd(xv3rij2, dx2ij);
					tmp2 = _mm_sub_pd(dv2ij, tmp2);
					__m128d j2ij = _mm_mul_pd(mimjrij3, tmp2);
					j2i =_mm_sub_pd(j2i, j2ij);
					#endif
				}

				#ifndef ACC_OFF
				_mm_store_pd(tmp, f0i);
				f0[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, f1i);
				f1[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, f2i);
				f2[i]+=tmp[0]+tmp[1];
				#endif
				
				
				#ifndef JERK_OFF
				_mm_store_pd(tmp, j0i);
				j0[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, j1i);
				j1[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, j2i);
				j2[i]+=tmp[0]+tmp[1];
				#endif
				
				#ifndef POT_OFF
				_mm_store_pd(tmp, poti);
				// i==j problem :: tempory :: substitute when will separate potential....
				// NB! do not use eps2=0 !, use small eps2 like 1e-10 or so 
				if (i==_j11) tmp[0]+=2*SQR(m[i])/sqrt(eps2);	
				
				pot[i]+=tmp[0]+tmp[1];
				#endif
			}
		}
		
	}	

	
	pthread_barrier_t barrier;	
	
	struct TH {
		FJPF *ff;
		double *f0,*f1,*f2, *j0,*j1,*j2, *pot;
		char *f0p,*f1p,*f2p, *j0p,*j1p,*j2p, *potp;
		pthread_t pth;
		uint id;
		uint i1_, j1_, i2_, j2_;	
		uint n1, n2;
	}*th;
	
	FJPF(){
		n_core=0;
		N=n=0;
	}
		
	uint get_i02(lluint nij, uint n){
		//if (nij<1) printf("get_i()::error nij<1\n");
		llreal D = (llreal)((lluint)n*(n-1) - 2*nij) + 0.25;
		llreal k = (0.5*(2*n-1)-sqrt(D));
		return fabs(ceil(k-1));
	} 
			
	uint get_j02(uint nij, uint n, uint i){
		uint k = i+1;
		uint nijn = ((lluint)k*(k-1))/2 + (lluint)k*(n-k);
		return n-1-(nijn-nij);
	}

	uint get_i1(uint tnij_max, uint thid, uint n){
		return get_i02(tnij_max*thid+1, n);
	}

	uint get_j1(uint tnij_max, uint thid, uint n, uint i){
		return get_j02(tnij_max*thid+1, n, i);
	}
			
	static void set_ncoreopt_nijmaxmin(luint nij, uint &n_core, lluint &nij_max, lluint &nij_min){
			
		//lluint nij = ((lluint)n*(n-1))/2;
			
		lluint n_cc_min = nij/n_core;
		if (nij%n_core!=0) n_cc_min++;

		uint n_core_opt = nij/n_cc_min;		
		if (nij%n_cc_min!=0) n_core_opt++;
			
		n_core = n_core_opt;
		nij_max = n_cc_min;
		n_core>1 ? nij_min = nij - n_cc_min*(n_core_opt-1) : nij_min = 0;
	
	}
	
	// set work coordinates
	void set_work_cor(uint thid, uint n_core, uint n, uint N, lluint tnij_max, uint &i1, uint &j1, uint &i2, uint &j2){
			
		if (thid==0){
			i1 = 0;
			j1 = 1;
		} else {
			i1=get_i02(tnij_max*thid+1, N);
			j1=get_j02(tnij_max*thid+1, N, i1);
		}
		
		if (thid<n_core-1){
			i2=get_i02(tnij_max*(thid+1), N);
			j2=get_j02(tnij_max*(thid+1), N, i2);		
		} else {
			i2 = n-1; // or n-2 ?
			j2 = N-1;
		}

	}
		
	// corrects work coordinates: we need j1 to be even
	void correct_work_cor(uint N, uint &i1, uint &j1, uint &i2, uint &j2) {	// makes j1 even 
			
		if ((j2!=N-1)&& (j2%2==0)) j2++;
				
		if (j1%2!=0) 
			if (j1!=N-1) {
				if (j1!=i1+1) j1++;
			} else {
				j1=i1+2;
				i1++;
			}
	}
	
	// set sum coordinates
	void set_sum_cor(uint id, uint n_core, uint n, uint &i1, uint &i2) {
		
		if (n<n_core) n_core = n; //  reduce paralellization if n<n_core 
		
		if (id>=n_core) { i1=1; i2=0; return; } // do nothing if n>n_core and id > n
	
		uint di = ceil((double)n/n_core); // max work length
		
		i1 = id*di;
		
		if (i1>=n) { i1=1; i2=0; return; } 
			
		n-i1>di ? i2=i1+di-1 : i2 = n-1; // please think by yourself :)
		
		
	}
	
			
	static void* thread_luncher(void* any){
		((TH*)any)->ff->thread(((TH*)any)->id);
	}
	
	virtual void thread(uint id){
		
		TH thi = th[id];
				
		set_work_cor(id, n_core, n, N, nij_max, thi.i1_,thi.j1_,thi.i2_,thi.j2_);
		
		correct_work_cor(N, thi.i1_,thi.j1_,thi.i2_,thi.j2_);
		
		ff(n, N, 
				thi.i1_,thi.j1_,thi.i2_,thi.j2_, 
				m, 
				x0,x1,x2,
				v0,v1,v2,
				thi.f0,thi.f1,thi.f2,
				thi.j0,thi.j1,thi.j2,
				thi.pot, eps2);
		
		
		pthread_barrier_wait(&barrier);
		
		// combine forces, 
		set_sum_cor(id, n_core, n, thi.n1, thi.n2);
		
		double f0i,f1i,f2i, j0i,j1i,j2i, poti;
		
		for (uint i=thi.n1; i<=thi.n2; i++)	{
			
			f0i=f1i=f2i=j0i=j1i=j2i=poti=0;	
			for (uint j=0; j<n_core; j++) {	
				
				#ifndef ACC_OFF
				f0i+=th[j].f0[i]; 
				f1i+=th[j].f1[i];	
				f2i+=th[j].f2[i];
				#endif
				
				#ifndef JERK_OFF
				j0i+=th[j].j0[i]; 
				j1i+=th[j].j1[i];	
				j2i+=th[j].j2[i];
				#endif
				
				#ifndef POT_OFF
				poti+=th[j].pot[i]; 
				#endif
			}
			
			#ifndef ACC_OFF
			f0[i]=f0i;
			f1[i]=f1i;
			f2[i]=f2i;
			#endif
			
			#ifndef JERK_OFF
			j0[i]=j0i;
			j1[i]=j1i;
			j2[i]=j2i;
			#endif
			
			#ifndef POT_OFF
			pot[i]=poti; 
			#endif
		}		
		
	}
	
	virtual void ff(uint n, uint N, double *m, 
						double*x0,double*x1,double*x2,
						double*v0,double*v1,double*v2,  
						double*f0,double*f1,double*f2, 
						double*j0,double*j1,double*j2, 
						double*pot, double eps2, uint n_core) {

		tr.set();
		
		if (this->N!=N || this->n_core!=n_core)   // recreate x0s.., th[i].f0.. if N or n_core changed
			if (this->n_core)
				for (uint i=0; i<this->n_core; i++) {
					
					#ifndef ACC_OFF
					delete [] th[i].f0p;
					delete [] th[i].f1p;
					delete [] th[i].f2p;
					#endif
					
					#ifndef JERK_OFF
					delete [] th[i].j0p;
					delete [] th[i].j1p;
					delete [] th[i].j2p; 
					#endif
					
					#ifndef POT_OFF
					delete [] th[i].potp; // delete old array
					#endif
				}

		if (this->n_core!=n_core)	{
			if (this->n_core) delete [] th;
			th = new TH[n_core];
			
			for (uint i=0; i<n_core; i++) {
				th[i].id = i;
				th[i].ff = this;
			}
			pthread_barrier_init(&barrier, NULL, n_core);
		}

		if (this->N!=N || this->n_core!=n_core) 
			for (uint i=0; i<n_core; i++) {
			
				for (uint i=0; i<=4; i++) {
					char *testp; 
					double *test  = (double*)new_align(N*14, 8, 16, testp);
					delete [] testp;
				}
				#ifndef ACC_OFF
				th[i].f0 = (double*)new_align(N+2, 8, 16, th[i].f0p); //*10 -performance reason
				th[i].f1 = (double*)new_align(N+2, 8, 16, th[i].f1p);
				th[i].f2 = (double*)new_align(N+2, 8, 16, th[i].f2p);
				#endif
				
				#ifndef JERK_OFF
				th[i].j0 = (double*)new_align(N+2, 8, 16, th[i].j0p); //*10 -performance reason
				th[i].j1 = (double*)new_align(N+2, 8, 16, th[i].j1p);
				th[i].j2 = (double*)new_align(N+2, 8, 16, th[i].j2p);
				#endif
				
				#ifndef POT_OFF
				th[i].pot = (double*)new_align(N+2, 8, 16, th[i].potp);
				#endif
			}
		
		
		this->n_core = n_core;
 		this->N = N;
		this->n = n;
		this->m = m;
		
		this->x0 = x0;
		this->x1 = x1;
		this->x2 = x2;
		
		this->v0 = v0;
		this->v1 = v1;
		this->v2 = v2;
		
		this->f0 = f0;
		this->f1 = f1;
		this->f2 = f2;
		
		this->j0 = j0;
		this->j1 = j1;
		this->j2 = j2;
		
		this->pot = pot;
		
		this->eps2 = eps2;
		
		nij = (luint)N*n-((lluint)n*(n+1))/2;
		set_ncoreopt_nijmaxmin(nij, n_core, nij_max, nij_min);
		
		for (uint i=0; i<n_core; i++)
			pthread_create(&th[i].pth, NULL, &thread_luncher, (void*)&th[i]);
	
		for (uint i=0; i<n_core; i++)
			pthread_join(th[i].pth, NULL);
		
		tr.get();
	}

	void print(){
				
		uint op_nnij=13;
		uint op_nn2ij=13;
		
		#ifndef ACC_OFF
		op_nnij+=6;
		op_nn2ij+=9;
		#endif
		
		#ifndef POT_OFF
		op_nnij+=2;
		op_nn2ij+=3;
		#endif
		
		#ifndef JERK_OFF
		op_nnij+=22;
		op_nn2ij+=25;
		#endif
		
		
		lluint nn2ij =((lluint)n*(n-1))/2;
		lluint nnij = (luint)N*n-(luint)n*n;
	
		double op1 = (double)((op_nnij)*nnij+(op_nn2ij)*nn2ij)*1e-9;
		double op2 = (double)((op_nnij+18)*nnij+(op_nn2ij+18)*nn2ij)*1e-9; // 18 for sqrt
	
		printf("\n");
		printf("\033[22;34mop_nnij:\033[22;30m %u \033[22;37m\033[22;30m\n", op_nnij);
		printf("\033[22;34mop_nn2ij:\033[22;30m %u \033[22;37m\033[22;30m\n", op_nn2ij);
		
		printf("\033[22;34mN:\033[22;30m %u \033[22;37m\033[22;30m\n", N);
		printf("\033[22;34mn:\033[22;30m %u \033[22;37m\033[22;30m\n", n);
		printf("\033[22;34mncore:\033[22;30m %u \033[22;37m\033[22;30m\n", n_core);
		printf("\033[22;34mTime:\033[22;30m %g \033[22;37msec\033[22;30m\n", tr.sin);
		
		printf("\033[22;34mPerf1:\033[22;30m %g \033[22;37mGflops // (19+2+22)*nnij+(22+3+25)*nn2ij\033[22;30m\n", op1/tr.sin);
		printf("\033[22;34mPerf2:\033[22;30m %g \033[22;37mGflops // (18+19+2+22)*nnij+(21+19+3+25)*nn2ij 19-sqrt\033[22;30m\n", op2/tr.sin);
		
		printf("\n");
	}
	
	void print_P(uint nprint){
		
		for (uint j=0; j<=nprint; j++)  {
		
			uint i = (N-1)*double(j)/nprint;
			
			printf("[\033[22;34m%05d\033[22;30m]:  ", i); 
			
			#ifndef ACC_OFF
			printf("% .15e % .15e % .15e  ", f0[i]/m[i], f1[i]/m[i], f2[i]/m[i]);
			#endif
		
			#ifndef JERK_OFF
			printf("% .15e % .15e % .15e  ", j0[i]/m[i], j1[i]/m[i], j2[i]/m[i]);
			#endif
			
			
		
 			#ifndef POT_OFF
			printf("% .15e ", pot[i]);
			#endif
			
			printf("\n");
		
		}
	}


};

}
#endif
